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ABSTRACT 

We present a powerful new tool for fitting broadband gamma-ray burst afterglow data, which can be used to 
determine the burst explosion parameters and the synchrotron radiation parameters. By making use of scale 
invariance between relativistic jets of different energies and different circumburst medium densities, and by 
capturing the output of high-resolution two-dimensional relativistic hydrodynamical (RHD) jet simulations in 
a concise summary, the jet dynamics are generated quickly. Our method calculates the full light curves and 
spectra using linear radiative transfer sufficiently fast to allow for a direct iterative fit of RHD simulations to 
the data. The fit properly accounts for jet features that so far have not been successfully modeled analytically, 
such as jet decollimation, inhomogeneity along the shock front and the transitory phase between the early time 
relativistic and late time non-relativistic outflow. As a first application of the model we simultaneously fit the 
radio. X-ray and optical data of GRB 990510. We not only find noticeable differences between our findings for 
the explosion and radiation parameters and those of earlier authors, but also an improved model fit when we 
include the observer angle in the data fit. The fit method will be made freely available on request and on-line 
at http : //cosmo . nyu . edu/af terglowlibrary In addition to data fitting, the software tools can 
also be used to quickly generate a light curve or spectrum for arbitrary observer position, jet and radiation 
parameters. 

Subject headings: gamma-rays: bursts - hydrodynamics - methods: numerical - relativity 



1. INTRODUCTION 

Gamma-ray bursts (GRBs) are short intense flashes of 
gamma radiation produced by cataclysmic ste llar ev ents such 
as the collapse o f the core of a massive star (IWooslevHT993l: 
iPaczvnskil [l998t iMacFadven & WooslevI \l99W or a neu- 
tron star-neu t ron star or neutron star- black hole merger (e.g. 
lEichleret all 119891: IPaczvnskil 1 1 99 Ih . During these events 
a coUimated relativistic outflow is produced which sweeps 
up the matter surrounding the GRB. Regardless of the orig- 
inal mass content or launching mechanism of this outflow 
(be it a fireball, Meszaros & Rees 1997, or a Poynting flux 
dominated jet, e.g. .Drenkhaha .2002). the expanding blast 
wave will sweep up circumburst matter and will eventually 
start to decelerate. As the blast wave shocks the circum- 
burst medium, broadband synchrotron radiation is produced 
by shock-accelerated electrons, giving rise to an afterglow 
signal that can be observed for up to days at X-ray and op- 
tical frequencies, and for up to years at radio frequencies. 

Three kinds of parameters determine the shape of the ob- 
served afterglow light curves. First, the shock dynamics are 
set by the explosion energy and circumburst density. A second 
set of parameters captures the physics of synchrotron emis- 
sion from shock-accelerated electrons. Finally, the observed 
flux depends on the parameters defined by a given observa- 
tion: frequency, time and observer angle. 

Ever since the first afterglows were discovered (iCosta etaTI 
[1997; Groot et al. 1997), models based on synchrotron radi- 
ation from a decelerating relativistic blast wav e have been 
successful in describing the broadband data (e .g. IWiier s et al. | 



T997 [lSari et aT][T998l:lWiiers & Galamal[T999 : Granot & S^ 



20021) . The synchrotron spectrum is typically described as 



consisting of a number of connected power law regimes, with 
the critical frequencies connecting the regimes shifting over 



time and being determined by the basic spectral shape of syn- 
chrotron emission, synchrotron self-absorption and cooling of 
the shock-accelerated electrons. In order to accurately model 
the late time afterglow emission in the radio, afterglow mod- 
els based on a non-relativistic rathe r than rela tivistic blast 
wave have been applied as well (e.g. iFrail et aLl lToOO). Both 
the early time ultra-relativistic and late time non-relativistic 
stages of the evo lution of the shock are self- similar, with so- 
lutions given bv jBlandford & McKe ^ (1976'), her eafter BM, 
and Sedov ( 195^, von Neumann ^on Neumamij (119611) and 
ITayloil fll950 ). hereafter ST, respectively. The intermediate 
stage can be approximated (e.g. Rhoads 1999; Huang et al] 
[1992), but this stage is complicated by the dynamics of jet 
decollimation. When jet decollimation is taken into account, 
a homogeneous jet surface that w idens with t he comoving 
speed of sound is often assumed (Rhoads"1999^, while more 
recent jet spreading models (Granot & Piran 2011 ) do not take 
the radial structure of the outflow into account. The jet nature 
of the outflow ultimately reveals itself as a break in the light 
curve, the jet break, and a subsequent steepening of the light 
curve slope. The ph ysics of afterglow jets has been r eviewed 
in e.g. Piran (20()l; lMeszarosl(l2006l) : lGranoti(l2007l) . 

Purely analytical models are severely limited in that they 
do not accurately capture many features of the jet dynam- 
ics (such as the aforementioned jet spreading and decelera- 
tion) and radiation. The simplifications inherent in a purely 
analytical approach lead to diverging predictions for a range 
of features such as the observed shape of the jet break, the 
size and shape of the counterjet (which was launched away 
from the observer and is only seen at late times, when rel- 
ativistic beaming of the emitted radiation plays a lessened 
role), the nature and duration of the transition to the non- 
relativistic phase and the effect of the orientation of the jet 
with respect to the observer. To gain a better understanding of 
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these aspects, various authors have performed numerical rel- 
ativistic hydrodyn amics (RHP) simulations of afterglow jets , 
in one dimension fKobayashi et al."1999';'Downes et al."2002t 

fimica et al. 2009; Van Eerten et al. 2010b), two dimensions 
Iranot et al. 2001; Zhang & MacFadyen 2006; MeHani et aD 
2007; Ramirez-Ruiz & MacFadyen 2010; Van E erten et al] 
2011; Wygoda et al.i [201 lb and occasionally even three 
ICannizzo et aL, 2004 1). Over the past decade, these simula- 
tions have steadily increased in accuracy, mainly through the 
use of adaptive-mesh refinement (AMR) techniques, which 
locally increase resolution during simulation where needed, 
and are capable of resolving the six orders of magnitude dif- 
ference between the initial width of the thin relativistic shell 
and the late time outer radius of the decelerated jet. Re- 
cent high-resolution simulations have shown that relativis- 
tic jets spread sideways significantly slower than analyti- 
cally expected and have a strongly inhomogeneous shock 
front (Zhang & M acFad yen 2009; Meliani & Keppens 2010; 
Van Eerten et al .l 1201 ll Van Eerten & MacFadyen 201 la; 
Wygoda et al...201 lb . Furthermore, they show that the tran- 
sition from rel ativistic to non-relativ i stic expansion is a very 
slow process (" Zhang & MacF adyen' 20 091 : IVan Eerten et all 
HoiOb); and that jet orientation stro ngly affects the Jet break 
even for small observer angles (I Van Eerten et al.l l2010ah . 
while for observers both on and off axis the observed jet-break 
time differs bet ween different frequenc ies due to synchrotron 
self-absorption dVan Eerten et al.i r201 1*). 

A RHD jet simulation can be combined with a numerical 
synchrotron radiation calculation to yield a powerful tool to 
predict the evolution of the observable broadband afterglow 
spectrum in detail. The weaknesses of simplified analytical 
models are thus avoided and local changes in fluid structure 
and arrival time effects are correctly accounted for This cal- 
culation can be performed in different ways, for example by 
summing over the emitted power of all fluid cells in the sim- 
ulation in the case of an op t ically thin fluid (Downes et al] 
I200I iNakar & Grano3l200l IZhang & MacFadyen. 2009.) or 
by fully solving the Unear radiativ e transfer equations in- 
cluding synchrotron self-absorgtion dVan Eerten et aljr2010bl: 
IVan Eerten & MacFadyen 201 iB- 

An obvious drawback of simulation based light curves 
compared to analytically calculated light curves is that cal- 
culating the former is a time consuming process. A full 
jet simulation takes several thousand CPU-hours to com- 
plete. A purely analytical light curve, on the other hand, 
can be calculated almost instantaneously and can therefore 
be applied to iterative model fitting, where the procedure 
of minimizing requires at least thousands of light curve 
calculations with slightly differing explosion and radiation 
parameters. In this paper, we present a new method to 
use simulation results directly as a basis for iterative fit- 
ting of broadband data, which closes the gap between sim- 
ulations and analytical models. This method should prove 
useful for further constraining the physics of gamma-ray 
burst afterglows, thereby obtaining clues about the nature 
of the progenitor and the burst environment. This provides 
more accurate predictions for future surveys inclu ding LO- 
FAR (Rottgering 2006), SKA dCarilli & RawHngs 2004) or 
ALMA (Wootten 2 003i) , an d indirectly benefits gravitational 
wave predictions for L IGO (lAbbott et al.ll2009h and VIRGO 
(lAcernese et al.ll2008h . where GRBs are potentially observ- 
able as electro-magnetic counte rparts (Nakar & Piran 2011; 
I Van Eerten & MacFadvenir201 Ibl) . Finally, our method helps 



to establish a baseline for studies of t he effect of more detailed 
models of the m icrophysics (see e.g. IVan Eerten et aLll2010bl ; 
iPanaitescu et alJl2006l:lFrigas et al.ll20l'ilK 

This paper is structured as follows. First we briefly describe 
the numerical settings and code used for our RHD jet sim- 
ulations of jets expanding into a homogeneous circumburst 
medium in ^ Our approach is made possible by two proper- 
ties of decollimating and decelerating relativistic jets starting 
from a BM solution. First, the jet evolution is scale-invariant 
under rescaling of both explosion energy and circumburst 
density, which we discuss in ^ Secondly, for a given initial 
opening angle, the 2D fluid profile of the blast wave evolves 
smoothly from a relativistic and purely radial outflow to a 
non-relativistic and spherical outflow. This implies that both 
the radial and lateral structure of the flow can be captured with 
sufficient accuracy by a low resolution grid with specialized 
coordinates that can be determined a posteriori once the ra- 
dial and angular extent of the jet at each moment in lab frame 
time are known from the high-resolution simulation. The dy- 
namical evolution and condensed low-resolution description 
of jets with various opening angles are discussed in ^ In 
^ we describe how the simulation results have been imple- 
mented in a broadband fitting code, which we apply to a case 
study, i.e. GRB 990510, in g6l We discuss our results in ^ 

The source code of the broadband fit code will be 
made freely available on request or for download from 
http : / / cosmo . nyu . edu/ afterglowlibrary. It 
can be run both on a single core and in parallel, and allows 
the user to either quickly generate light curves and spectra for 
arbitrary explosion parameters, or, when provided with a data 
set of observed fluxes in mJy, to perform a full broadband fit. 

2. NUMERICAL JET SIMULATIONS 

For this study a total of 19 jet simulations in 2D 
have been performed using th e Relativistic Adaptiv e Mesh 
(RAM) parallel RHD code ( IZhang & MacFadvenI ^006). 
The code employs the fifth-order weighted essentially non- 
oscillatory (WENO) scheme (Jiang & Shu 1996) and uses 
the PARAMES H AMR tools dMacNeice et all 12000) from 
FLASH 2.3 (.Fryxell et al.ll2000l) . For aU jet simulations the 
BM solution for an adiabatic impulsive explosion is used in 
spherical coordinates to set the initial conditions. Instead of 
the full spherical solution, a conic section is used that is trun- 
cated at a different fixed opening angle for each simulation. 
The opening angles are listed in table [U along with the jet 
energy Ej for each simulation. 

The jet energy Ej (the total for both jets) relates to the 
isotropic equivalent energy Eiso according to 

Ej = £■,■„,(! -cos 6*0) ~ EisoOl/l. (1) 

All jets expand into a homogeneous medium with number 
density «o = 1 cm"-* (mass density po = 1 x nip g cm"-', in 
terms of the proton mass nip), and have an isotropic equiv- 
alent explosion energy = 6.25 x 10^' erg; note that due to 
the scale invariance of the simulations with respect to no and 
Eiso these can be scaled afterwards to represent arbitrary val- 
ues (see Section[3]l. All jets start at time ti, with fluid Lorentz 
factor = 25 directly behind the shock, ensuring that for all 
simulations 7^, > I/^q. At this point the edges of the jet have 
not yet come into causal contact and lateral spreading has not 
yet set in, which allows us to use the spherically symmetric 
BM solution as the starting point. 
The initial outer radii of the jets are given by Rb= 1 .3 102 x 
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TABLE 1 

Opening angles in radians and jet energies Ej in erg for 
each simulation. 





6*0 (rad) 


Ej (erg) 


1 


0.045 


6.328 




2 


0.05 


7.812 


IQ*» 


3 


0.075 


1.758 


10*" 


4 


0.1 


3.125 


1049 


5 


0.125 


4.883 


1049 


6 


0.15 


7.031 


1049 


7 


0.175 


9.570 


1049 


8 


0.2 


1.250 




9 


0.225 


1.582 


1050 


10 


0.25 


1.953 


lO^o 


11 


0.275 


2.363 


1050 


12 


0.3 


2.813 


1050 


13 


0.325 


3.301 


1050 


14 


0.35 


3.828 


1050 


15 


0.375 


4.395 


1050 


16 


0.4 


5.000 


1050 


17 


0.425 


5.645 


1050 


18 


0.45 


6.328 


1050 


19 


0.5 


7.812 


1050 



lO'^ cm, determined from th and 7/, through Rh = cf/,(l - 
1/167^) (eq. 26 from BM, with c the speed of light). The 
initial time is determined from 7^,, and no using 

llEiso=l6TTmpnQ^lc^tl, (2) 

which expresses conservation of total energy (eq. 43 from 
BM). 

For all simulations the stopping time is determined accord- 
ing to 

/ E V 
f/=10xfNR = 9700fy^j days (3) 

The time Jnr marks the time when the jet is analytically ex- 
pected to transition from relativistic to non-relativistic flow, 
and is determined from comparing the initial explosion en- 
ergy to the total rest mass energy of the swept-up matter 
(jPiran 2005). Numerical simulations have shown that in prac- 
tice this transition takes more time to complete, which is 
why we have chosen to continue our simulations until two 
times the transition duration 5 x Tnr found numerically by 
[Zhang & MacFadyen (2009). For Ei^o = 6.25 x 10^' erg and 
7Z, = 25, it follows for all 19 simulations that tj, = 4.37 • lO*" 
s = 50.6 days and tf = 3.33 • 10** s = 3849 days; again, these 
values change under rescaling of Ejso or no- 

All simulations use an equation of state with the adiabatic 
index as a function of comoving density and pressure chang- 
ing smoothly fro m 4/3 for a relativistic fluid to 5/3 for a non- 
relativ istic fluid jZhang & MacFadveiil 120091; iMignone et alj 
l2005h . 

2.1. Resolution and Refinement 
The initial widfli of flie shell is ARb ^i/27^ « 1.05 x 
10'"^ cm. In order to correctly resolve this initial width, 
we have set the initial peak refinement level to 15, which 
for a grid running from 0.01 x 7?^ to c x tf and 384 radial 
cells at the base level, implies a smallest radial cell size of 
Sr = 1.58 X 10'^ cm, since each increase in refinement level 
double the effective resolution. There are 32 base level cells 
in the angular direction, so that 69 = 9.6 x 10"^ rad. Over 



time, the peak refinement level is gradually decreased un- 
til a peak level of 9 (the sa me approach has been applied in 
IZhang & MacFadvenil2009l) . 

In addition to the global change in peak refinement level, 
a number of additional manual derefinement strategies have 
been employed in order to prevent the simulation from de- 
voting too much of its calculation time and memory to re- 
solving the boundary and non-relativistic sideways shock be- 
tween the interstellar medium (ISM) and empty region far 
behind the shock front, as well as the Kelvin-Helmholtz in- 
stabilities that arise in the inner low density regions of the 
shock due to velocity shear at the edge. These regions have 
no relevance for the jet dynamics or the observed radiation, 
which is produced close to the front of the shock. The shock 
front and its sideways expansion are fully resolved. The ad- 
ditional manual derefinement settings are an increasing inner 
radius behind which derefinement is enforced. This region 
expands as r = 1.2 x 10'^(f/fi) until 7 = 2 (according to BM) 
and then stops. Regions where the fluid number density is 
below 0.75 no have their peak refinement level reduced by 6. 
The peak level is also reduced by 6 where the local 7 < 1.5, 
but is never allowed to drop below 9 for derefinement based 
on this criterion. Finally, the peak level is reduced by 4 where 
the local 7 < 3, but is never allowed to drop below 1 1 for 
derefinement based on this criterion. 

3. SCALE INVARIANCE OF THE JET 

What is straightforward, but perhaps obscured by the rich- 
ness of features in the resulting light curves, and what has so 
far not been utilized in afterglow modeling, is that the full 
2D evolution of the jet is invariant under scaling of the initial 
explosion energy and circumburst medium density: indepen- 
dent of self-similarity, a more energetic jet (or a jet in a less 
dense environment) goes through exactly the same evolution- 
ary stages as a jet with lower energy (or higher circumburst 
density), albeit that each stage occurs at a later time and larger 
radius. 

The scalings can be understood as follows. The initial set- 
up of the problem is determined completely by a limited num- 
ber of parameters: £,50, po (= no x nip, with nip the proton 
mass), 9o and c. The initial Lorentz factor 7^ is not included 
in the list because its precise value is arbitrary as long as 
7ft > l/6'o. Only a limited number of independent dimension- 
less combinations of these parameters is possible, such as: 

A=-, B=^, 9, 9,. (4) 
ct pqv^ 

Any dimensionless quantity that describes the local fluid con- 
ditions or global evolution of the jet, such as Lorentz fac- 
tor 7, density p/po or 6*95 (the angle with respect to the 
jet axis within which 95% of the jet energy is contained) 
can be expressed as a function of these parameters, i.e. 
^{r,t,9,Eiso,po,9o) j{A,B,9,9()). Note that at both very 
early times (no lateral flow) and at very late times (spheri- 
cal flow) there is no dependency on 9 or 9o, which together 
with the fact that A becomes a constant value both in the ultra- 
relativistic BM and non-relativistic ST limits (i.e. A — > 1 and 
A 0, respectively), accounts for the self-similarity of these 
limiting cases, with j{r,t,9) jiB), etc. The parameter B 
is identical to the ST self-similarity variable and is cen- 
tral to the self-similarity of the BM solution via equation [2] 
above (where r ^ ct), which fixes 7^ and therefore the BM 
self-similarity variable x(B,j^(B)). 
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G. 0x10° 4.0x10= 2.0x10' 5.0x10= I.OxiO" 1.5x10° 

Log 10 tau 




6.0x10= 4.0x10° 2.0x10' S.OxlO'^ 1.0x10° 1.5x10° 

Fig. 1 . — Direct comparison between comoving number density n in cm"^ 
(top) and lab frame energy density r in units of nipC^ (bottom) profiles for 
jet simulations with 0q = 0.2 rad, no = 1 cm' "3 , and Ei,o = 5x 10^' erg (left) 
or £■,.„ = 5x10"' erg (right), drawn from lVan Eerten & MacFadvenI <2011b ). 
The snapshot times differ by a factor 100'/^ and for both snapshots the fluid 
Lorentz factor directly behind the shock is ~ 3.3. 

Even if we do not limit ourselves to either of the self-similar 
extremes and leave A and 9 in place, we have a set of coordi- 
nates A, B, 6 that is invariant under the following rescalings: 

Po = Vo, 

t'=(K/xy/^t. (5) 

The full scale-invariant hydrodynamics equations in terms of 
A, B, 9 are provided for completeness in appendix |B] The 
direct implication of this is that we can determine the value 
of any dimensionless quantity for an explosion with E'^^^ and 
P() simply by probing a simulation with parameters Eiso,p at 
t,r, rather than at t\r'. Quantities that are not dimensionless, 
such as density p and internal energy density e follow from 

P/Po = P'/P'o, e/ poc^ = e'/p'oc^, etc. 

Examples of scalings between jets with different Ejso and 
po are given in Figures [T] and |2] The comoving fluid density 
and energy density profiles are drawn from the simulations 




I.QxID^ 5,0x10* 5,0x10^ 1.0x10^ 1.5x10^ 
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Fig. 2. — Direct comparison between comoving number density n in cm"^ 
(top) and lab frame energy density t in units of iripC^ (bottom) profiles for 
jet simulations with 9q = 0.2 rad, £,.,0 = 5 X 10*" erg, and no = I X 10^^ cnC^ 
(left) or no = 1 cm"^^ (right), drawn from Van Eerten & MacFadvenI i201 Ihll . 
Unlike in the case of the two simulations of Figure [T] there is now a scaling 
factor A = 10^ between fluid quantities p and r for the different fluids as well. 
For both snapshots the fluid Lorentz factor directly behind the shock is ~ 4. 

presented in I Van Eert en & Ma cFadyeiil (1201 Ibl) . since the 19 
simulations listed in table [T] were set up to be unrelated via 
scaling. In the case of Figure [T] two simulations are com- 
pared for which A = 1, whereas A = lO'' in the case of Figure 
|2l In Figure[T] we set the outer radius of the left and right pan- 
els equal to their respective ct, resulting in both images being 
complete mirror images of each other, which confirms that 
both simulations are numerically indistinguishable. Although 
analytically equivalent, the two simulation runs in Fig. |2]are 
no longer numerically identical, which leads to minor differ- 
ences in (de-)refinement in the inner regions. However, those 
inner regions do not contribute to the observed flux, and thus 
have no observable effect on the light curves. 

4. LATERAL SPREADING AND JET DECELERATION 

We have plotted a number of features of a subset of 
the 19 jet simulations in Figures [3] and ID The jets fol- 
lo w the same evolution a s thos e of earlier of studies, such 
as Zhana & MacFadven ('20091); IVan Eerten & M acFadvenI 
(2011 a); Wygoda et al. (2011). The top plot of Figure[3]shows 
the evolution of the on-axis outer radius of the blast wave for 
jets with different 9q in the lab frame of the explosion. All 
jet simulations start out relativistically with identical 7/,, tb 
and blast wave radius Rb- Jets with a wide opening angle do 
not have a long spreading phase and undergo a single smooth 
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lab frame time t (days) 




lab frame time t (days) 



Fig. 3. — Top: evolution of the blast wave radius R over lab frame time for 
jet with 6»o = 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05 rad (top to 
bottom in both plots). Dashed green lines indicate asymptotic BM and ST 
predictions, both for a spherical blast wave and for 9q = 0.05 rad. Bottom: 
evolution of blast wave velocity fi-y for the same opening angle jets. Early 
time dashed green line indicates BM prediction of /3r7r oc r^^^, late time 
indicates /3r7i- oc r^/^ . 

transition from ct to 1 .15(EjP / poY^^ (where the nu- 
merical factor 1.15 follows from the ST solution with adia- 
batic index 5/3). The figure shows that for a very wide jet 
with 6*0 = 0.5, the transition time is well approximated by the 
crossing point of the asymptotes for spherical outflow. How- 
ever, this does not carry over to smaller angle jets, and the 
meeting point of the asymptotes for Ej{6o = 0.05) severely 
underestimates the turnover point. The reason for this is that 
for narrower jets there is also an intermediate phase, where 
the jet decelerates due to lateral spreading. Although not as 
ab rupt as originallv predicted (Rhoads 1999. also discu ssed 
in I Van Eerten & MacFadvenll201 latlWygoda et alJl20Tlh and 
occurring throughout the transrelativistic phase of the jet, this 
leads to an extended period of jet deceleration in excess of 
the asymptotic jet deceleration of the ST phase and adds a 
second turnover point in the evolution of jet radius and jet 
velocity. As the bottom plot of Figure [3] indicates, this inter- 
mediate phase does not follow a simple power law. A full pa- 
rameterization of the intermediate phase lies beyond the scope 
of this work and will not be required for model fitting based 
on the simulation results. In addition, due to inhomogeneity 



along the shock front, deceleration will be different for out- 
flows along different angles. 

The early time behavior visible in the bottom plot of Fig- 
ure [3] where the peak Lorentz factor initially drops below 
its expected value in the BM regime but then moves back 
to the BM asymptote, is due to the resolution of the simu- 
lations. In the BM solution, the blast wave is extremely thin 
(ARh ^ Rb/2'^l) and we have not been able to achieve full 
convergence in 2D at earl y times (the issue is id entical to that 
illustrated in Figure 5 of Zhang & MacFadver7 '2009). By us- 
ing the integrated values of the fluid quantities across a single 
fluid cell rather than the values at its central coordinate, we 
have ensured that all energy of the BM solution is accounted 
for during the initialization of the simulation. For this reason, 
the drop is only temporary. We emphasize that this is not due 
to lateral spreading of the jet, for if that were the case the peak 
Lorentz factor would not have been able to recover 

Lateral spreading and the inhomogeneity of the shock front 
are illustrated in Figure ID The plots show the time evolu- 
tion of various characteristic angles 6gs , djs , 650 of the outflow, 
defined as the angles within which a fraction 0.95,0.75,0.50 
of the volume-integrated rest frame energy density t is con- 
tained. The green dashed lines indicate the analytically ex- 
pected onset of lateral spreading, when 7^1 /Oq. Because 
the plots show Ogs,... rather than the outer edge angle 9 edge, 
these lines have been shifted according to 9(t) — > \/O350edge 

(or ^/Q.lSOedge or ^/Q.SQOedge), This means that if the simu- 
lations had followed the analytical estimate, the turnovers for 
all angles and fractions would have occurred along the green 
lines. Because the jets become strongly inhomogeneous along 
the shock front, the turnover is delayed for characteristic an- 
gles bounding smaller energy fractions: the jet energy stays 
concentrated near the jet axis even after the edges have be- 
gun to expand. Because narrower jets have smaller values 
of Ej, they will decelerate earlier, which leads to the behav- 
ior shown in the plots where ^95,... for small jets cross that 
of large jets with the same £,50. Jets with the same Ej do 
not show this effect (the characteristic angle curves for wider 
jets would shift sufficiently far to the left in the plots if their 
energies were downscaled to match the Ej of the narrowest 
jet). At late times all curves tend to the same fixed fraction of 
7r/2 regardless of initial energy, i.e. the jet becomes spherical 
and homogeneous along the shock front, but only for high en- 
ergy fractions is this limit actually reached in our simulations. 
The early time turnover behavior confirms the validity of our 
choice of starting 7^ > 6*0 (note that for the wide jets 7/, ^ Oq). 

5. "BOX"-BASED BROADBAND AFTERGLOW 
FITTING 

Each of the 19 simulations generates 3000 snapshots, vary- 
ing in data size between ^ 350MB at early times and ^ 40MB 
at late times (when lower peak refinement levels are utilized). 
Therefore it is currently not practically possible to load the 
complete output for a single simulation into computer mem- 
ory at once, let alone the combined output of all 19 simula- 
tions. However, if we wish to use the simulations as a basis 
for iterative model fitting, we will need to be able to quickly 
access the fluid state at any requested point in time and space. 
We therefore need to summarize the simulation results in a 
way that adequately captures all aspects of the outflow but 
occupies only a relatively small amount of computer memory. 
In this study we have aimed for < 1GB in total, but the desired 
size in memory will be hardware-dependent. 



6 




10^ 10^ 
lab frame time (days) 




10^ 10^ 
lab frame time (days) 




10^ 10^ 
lab frame time (days) 



Fig. 4. — Lateral spreading of jet simulations with opening angles Bq = 
0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05 rad (top to bottom in each 
plot). The top plot shows the evolution of Sgs, defined as the angle within 
which 0.95 of the volume-integrated rest frame energy density t is contained. 
Center plot shows the same for 675 (0.75 of energy) and the bottom plot 
for 650 (0.50 of energy). Dashed green lines indicate when 7 = \/9q for 
the different opening angles (shifted to correct for fractions < 1 ) and should 
therefore be compared to the initial opening angles of the jets, not to 9xx 
angles at any given time. The horizontal dotted green lines indicate the values 
for the hmiting angles for the case of a homogeneous spherical blast wave. 



The next three subsections describe the further technical as- 
pects of producing Ught curves from the summarized 'box', 
rather than the 'grid' which originally contains the simula- 
tion data. In subsection 15 .41 we describe the details of the fit 
method. 

5.1. "Box " summary of simulations 

The phase space for each fluid variable is set by the vari- 
ables Eiso, «(), On, r, 9 and t and the local fluid state is fully 
specified by the fluid variables r, p, vv (radial fluid velocity), 
vg (angular fluid velocity), from which the other fluid quan- 
tities such as p (pressure), e or 7 can be derived. By stor- 
ing e explicitly along with p, v^, vq, while still including r, 
even though the latter is not needed for the radiative trans- 
fer calculation, we ensure that we can explore the full fluid 
state at each point in time and space without having to use 
the equation of state. This implies that we will have to store 
five fluid variables along with the coordinates and sizes of the 
fluid cells they refer to (i.e. 9 quantities in total). The scale 
invariance discussed in section [3] implies that we only need 
to store a single value for Eiso and no- As mentioned already, 
we have used 19 possible values for 6q. We calculate light 
curves for intermediate values of by interpolating the ap- 
propriate resu lts f rom the 19 simulations, as we demonstrate 
in subsection 15.31 For the r, coordinates we use 100 en- 
tries each. For t we also use 100 entries rather than 3000, 
and as w ith the small subset of Oq, we demonstrate in subsec- 
tion 15.31 that this number is sufficient and that the full light 
curve can be calculated using interpolation. All together, we 
now have the following number of entries in our summary 
of the fluid evolution: Ei^,, x x 9o x r x 9 x t x variables = 
1 X 1 X 19 X 100 X 100 X 100 X 9 = 171,000,000. We will 
refer to this summary as the "box", in order to distinguish it 
from grid based fluid profiles from the simulations. 

The reason that we only need as few as 100 cells in the r 
and 9 directions has already been demonstrated partially by 
Figures [3] and |4] in the preceding section. These Figures illus- 
trate that, although high-resolution simulations were required 
in order to accurately determine the blast wave radius and lat- 
eral extent at each point in time, the evolution from BM to 
ST profile itself is smooth. This means that, once the large- 
scale properties of the outflow are known, it is possible to use 
this knowledge to select at each point in time new coordinates 
such that the key features of the outflow are resolved while 
ignoring parts of the outflow. 

The box 9 at each point in time will not cover the entire grid 
but only runs from to Omax = 6*99 /0.99 (i.e. slightly over the 
the outer angle within which 99 % of the volume integrated r 
is contained). Anything outside these angles is either ISM or 
contributes a negligible amount of radiation. The lateral cells 
within this region are evenly spread. Alternatively we could 
have set the cell coordinates using 6*00, ^'oi, etc., but there is no 
significant difference in the resulting light curves. The radial 
domain runs from to the outer limit Rmax of the blast wave 
at each box value of 9. Because the unshocked ISM is at rest, 
the outer boundary of the shock wave is readily determined 
using the criterion that 7 > 1.000001. Even for extremely 
high resolution, this point will not exactly coincide with the 
peak of the blast wave. We therefore devote 10 cells to the 
region between R determined by the r peak of the blast wave 
and the first shocked ISM cell. Of the remaining 90 radial 
box cells, 80 are devoted to resolving the blast wave. Since 
we know from the BM solution that the blast wave width is 
Ai? w /?/127^ initially and from the ST solution that AR(xR 
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eventually, we analytically determine the width of the blast 
wave by AR = R/ll-y'^, with 7^ ee ^b(t /tt)~^l'^+ 1 (where the 
' 1 ' has been added relative to BM eq. 24 to obtain the correct 
asymptotic behavior). Note that this is only approximately 
correct on-axis due to 2D spreading and less accurate for the 
radial profile along the outer angles. This does not matter, 
however, as long as the approximation is sufficient to resolve 
the sharp feature of the blast wave. The final ten radial cells 
are spaced between the origin and the back of the shock. All 
radial cells are equally spaced within their respective region. 
Figure |5] shows radial fluid profiles for the 6*0 = 0.05 simula- 
tion for different times and angles. Profiles for other values of 
6q are similar 
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Fig. 5. — Comoving number densities n for tlie 6*0 = 0.05 rad simulation at 
lab frame times t = 57.6 days (top) and f = 186 days (bottom). At these times 
the BM solution predicts respectively 7 = 20.6 and 7 = 3.6 behind the shock 
front (see also Figure[3j. The fluid profiles from the grid at the indicated an- 
gles are drawn with lines, while the box data points are presented by circles. 
The angles represent values close to the jet center, edge and, in the bottom 
plot, halfway between center and edge. The grid values are taken exactly at 
the listed angles, while the box values are averages centered on these angles 
with 59 = 5.75 X 10^ rad and 69 = 4.4 X 10"^^ rad at early and late time re- 
spectively. This accounts for small differences between box and grid values. 
The large 9 profile at the left in the bottom plot differs from the others in that 
many box cells are used to resolve the steep drop in front of the blast wave 
as well, which illustrates that the peak values of n and r are not numerically 
identical in fluid simulations (in this case, the peak value of t occurs at larger 
radius than the peak value of n). 



5.2. "Box" interpolation 

The previous subsection refers to interpolation in t as well 
as in 6q. Even though 100 time snapshots are adequate to 
capture the dynamics of the jet, evaluating the linear radiative 
transfer integral at just these 100 emission times (that is, in 
addition to any evaluations of the BM solution at times before 
the initial simulation time) has been found to lead to noise 
in the light curve both at early times and during the rise of 
the counterjet at late times (corresponding to early emission 
times for the counterjet). The solution to this is to interpolate 
the fluid profiles between different emission times, something 
that is neither difficult nor time consuming when done based 
on the box snapshots. We have found for light curves with 6q 
values between those in Table [T]that interpolation at the fluid 
level is more effective than interpolation at the level of the 
light curve: calculating two light curves for adjacent tabulated 
6q values and interpolating between them will systematically 
shift the jet break time and post-break asymptote relative to 
their actual values for the intermediate 6q. Therefore, both the 
time and Oo interpolation occur at the same stage. 

In practice ^ 3000 interpolated times have been found to be 
more than sufficient to remove the numerical noise in the light 
curve. Implementing the Oo and t interpolations to obtain the 
local fluid state for requested fluid coordinates r,9,t works as 
follows: 

• There are four tabulated entries needed for the inter- 
polation, determined by the closest surrounding values 
around the requested 6q and t. For each entry, scale 
&MAX by first interpolating in f, then in 6q. 

• For all four entries, determine the appropriate Rmax{G), 
after the 6 coordinates' outer boundaries have been 
scaled to their new Omax- Interpolate this Rmax in t 
for both 9{) options separately. 

• Obtain the box fluid conditions at r/in/Xy/^, 9 for all 
four entries, applying both the 9max scaling and Rmax 
scalings. 

• Multiply all non-dimensionless quantities by A. 

• Interpolate the results first in f , then in 9o to obtain the 
final value for the fluid quantity. 

We use linear interpolations, which produces converged re- 
sults (see ^5.3l l. 

5.3. Light curves 

We calculate light curves and spectra from the box us- 
ing the same method that we ha ve used previously for 
grid-based light curve calculations (I Van Eerten et al.|[2010at 
Van Eerten & MacFadven 201l3). The dominant radiation 
mechanism is assumed to be synchrotron radiation and the 
broadband emission from each fluid cell is a given by a serie s 
of connected power laws similar to those in lSari et"al] (119981) . 
The linear radiative transfer equations are solved simultane- 
ously for a large number of rays. For each point in lab frame 
time f , the plane perpendicular to the direction of the observer 
and at a fixed distance Reds from the origin of the box (or 
grid), defined hy Reds / c = t-tohs, defines the area from which 
emission will arrive at exactly the same observer time tots- 
This plane is labeled th e equidistant surface (EDS, see also 
IVan Eerten et al.ll2010 b'). In earlier work we have employed 
a procedure analogous to AMR for dynamically changing the 
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number of rays through the EDSs that are followed simulta- 
neously. For the blast waves of this study all EDS refinement 
would have occurred near the center of the EDS (defined as 
the intersection of the line from the grid origin to the observer 
and the EDS), in order to resolve the early time BM profile 
and the refinement level would have gradually decreased out- 
wards. Since this is essentially equivalent to base 2 logarith- 
mic spacing, for the current study we use fixed logarithmic 
spacing between rays in the radial direction instead. The num- 
ber of rays is evenly spaced in the angular direction. This 
approach is somewhat faster than the dynamical refining and 
requires less memory for bookkeeping. In the final step the 
rays are integrated over to yield the observed flux. Flux, ob- 
served frequency and observer time are corrected for redshift 
z during the calculation. 

In Appendix IaI we give the exact expressions for the emis- 
sion and absorption coefficients that are calculated while solv- 
ing the radiative transfer equations through the evolving fluid. 
Their local values depend on a number of parameters that 
capture the microphysics behind the synchrotron radiative 
process as well as on the local fluid conditions. There are 
four such parameters: the power law slope p of the shock- 
accelerated electrons, the fraction of magnetic energy rela- 
tive to thermal energy, the fraction tg of downstream thermal 
energy density in the accelerated electrons, and the fraction 
of the downstream particle number density that participates in 
the shock-acceleration process. By performing broadband fits 
on afterglow data using the box-based fit method from this pa- 
per, these parameters can be determined from the fits for the 
first time using the full blast wave evolution. 

All simulations start at 7/, = 25. Before this time the out- 
flow can be described by a conic section of the spherically 
symmetric BM solution. Because the observed flux at a given 
observer time contains emission from a wide range of emis- 
sion times, initially including times for which 7 > 7/,, the BM 
solution is used directly to determine the initial emission and 
absorption coefficients. Probing the BM solution between 
7 = 200 and 7;, using 1600 logarithmically spaced emission 
times has been found to be more than sufficient to capture the 
early time emission. 

We have performed various tests to check the resulting light 
curves from box based calculations. First of all, we have com- 
pared box based light curves to simulation based light curves 
for the simulation underlying the box summary. An example 
is shown in Figure|6] At most, the difference between the two 
is on the order of a few percent (see bottom plots in figure). 
The exceptions are the early time light curves for observers 
far off-axis, when the box approach smoothens out numerical 
noise more than the direct simulation calculation does. 

In Figures |7] and [8] we show a comparison to results of 
earlier studies. The close match between the box radio and 
optical light curves (solid and dashed lines for on and off- 
axis observer respectively) and the grid based counterparts 
(thick dotted lines) in Figure |7]is remarkable given the differ- 
ences in the methods by which the y were obtained. The grid 
based light curves are drawn from IVan Eerten et al.l (l2010ah 
and are therefore based on a simulation with different refine- 
ment and resolution settings and strongly differing isotropic 
energy compared to the unsealed simulations of the current 
work. In addition to that, the earlier curves have been calcu- 
lated with a completely different radiation algorithm, where 
a summation was done over the emitted power of each fluid 
cell (which therefore excludes the possibility of synchrotron 
self-absorption), rather than by employing the linear radiative 
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Fig. 6. — Direct comparison between box light curves and simulation light 
curves for simulation 8 (60 = 0.2 rad) of this paper. In the top plot, box based 
curves are shown both for optical (solid lines) on-axis and off-axis (Bobs = 0.4 
rad) and radio (dashed lines) on-axis and off-axis. Off-axis light curves start 
lower and peak later than their on-axis counterparts. The simulation based 
curves are drawn as wide dotted lines. The bottom four plots show the relative 
differences between box and simulation for all four combinations, according 



transfer method used in the current work. On average (over 
the logarithmically spaced data points) the difference between 
the box light curves and the light curves from the earlier study 
shown in the figure is a factor 1.15, with the biggest difference 
(a factor 1.31) occurring for the optical light curves (both on 
and off-axis) around 400 days. The fact that the off-axis light 
curves at some point cross the on-axis light curves and tem- 
porarily show higher flux levels is a result of relativistic beam- 
ing: at its most extreme it leads to the prediction of orphan 
afterglows, where the on-axis light curve remains effectively 
invisible relative to the off-axis light curve for observers at 
very high angles. 

In Fig ure [8] w e show a compar ison to light curves from 
IVan Eerten & MacFa dyeiil (1201 Ibh . These were obtained 
from simulations with lower resolution compared to the cur- 
rent simulations, which started from 7/, = 10 rather than 7^ = 
25. This accounts for the early time differences up to ^ 20 
days. Early time differences aside, the average difference 
(over the logarithmically spaced data points) is a factor 1 .09, 
with the biggest difference (a factor 1.23) occurring at late 
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Fig . 7. — Direct comparison between box light curves and grid light curves 
from lVanEertenetalJ )2010ij . see Figure 1 of that paper). The simulation 
from that paper had Oq = 0.2 rad, Ej^o = lO'^ erg, no = I criT^, z = 0, di_ = 
10^^ cm, p = 2.5, Ee = O.I, eg = 0.1, = 1, and the resulting light curves 
are indicated by thick dotted lines. The solid curves show the box results 
for Bats = rad at lO' Hz (radio, top) and lO'* Hz (optical, bottom). The 
dashed curves show box results for 9„i,s = 0.4 rad, also at lO' Hz (top) and 
10*'' Hz (bottom). Because the radiation was c alculated by direc t summ ation 
of the emitted power of all fluid elements in IVan Eerten et alj J2010aD . no 
synchrotron self-absorption was included in that work and for the purpose 
of comparison we have therefore disabled synchrotron self-absorption in the 
box curves for this plot as well. 
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Fig . 8. — Direct comparison b etween box light curves and grid light curves 
from lVan Eerten & MacFadvenI 1,201 lb, see Figure 1 of that paper). We use 
case B from that paper with Oq = 0.4 rad, Eiso = 1.25 X 10"" erg, hq = 10"' 
cm"', z = 0,dL = 10^^ cm, p = 2.5, £<■ = 0.1, eg = 0.1, = 1, and the resulting 
light curves are indicated by thick dotted lines. The solid curves show the 
box results for 9oi„ = rad at 1.43 GHz (radio, top) and 4.56 X lO'** Hz 
(optical, bottom). The dashed curves show box results for Soj, = 0.8 rad, also 
at 1 .43 GHz (t op) and 4.56 X lO''* Hz (bottom). The emission coefficients in 
IVan Eerten & MacFadven 1 201 lb) were larger by a factor 3/2 and in order 
to allow for a direct comparison, the fluxes from lVan Eerten & MacFadvenI 
1201 Ibl) have been divided by 3/2. 

times for the two off-axis light curves. 

Finally, in Figure |9] we show a comparison between opti- 
cal and radio light curves based on a single opening angle 
simulation (0.2 rad) and light curves for the same opening 
angle reconstructed from interpolation between simulations 
with ^0 = 0.175 rad and = 0.225 rad. The figure shows that 




observer time (days) 

Fig. 9. — Direct compaiison between box-based light curves with = 0.2 
rad, where no interpolation with respect to jet opening angle has been applied 
(thick dotted curves, using only data from simulation 8 in Table [T) and box- 
based light curves at Oq = 0.2 rad created by interpolation from 9o = 0. 175 
and 00 = 0.225 rad (solid and dashed curves, using data from simulations 7 
and 9). 

interpolated light curves generally match the light curves cal- 
culated from 9() = 0.2 rad. In this particular example the great- 
est discrepancy occurs between the low radio curves around 
160 days, where the interpolated curve flux is briefly larger 
by a factor 1.2. Note however that the figure represents an ex- 
treme scenario that does not occur in practice: by artificially 
removing 6'o = 0.2 for the purpose of testing the interpolation, 
we have tested interpolation across A6'o = 0.025 rad, whereas 
in practice the largest possible difference S9o = 0.0125 rad. 
We conclude that the range of jet opening angles between 
9o = 0.045 and 6*0 = 0.5 rad is adequately represented by our 
sample of simulations plus interpolation. 

The numerical errors between different simulations and box 
based light curves given above should be compared to the 
difference between simulation/box based light curves on the 
one hand and analytically calculated light curves on the other 
hand. Although the latter light curves have no numerical noise 
or resolution errors by definition, they have systematic errors 
due to the simplifications in the underlying assumptions for 
dynamics and radiating region that are far lar ger than the nu- 
merical noise in the former In Van Eerten et a l. (2010a) simu- 
lation results are compared to different analytical models and 
differences up to an order of magnitude in flux level and in 
time (for specific features such as the off-axis moment of peak 
flux) are seen, especially in the transrelativistic phase. 

5.4. Fitting methods 

The method to quickly generate the observed flux at an arbi- 
trary observer frequency and time described in the preceding 
subsections allows for iterative fitting of the simulation based 
afterglow model of a decelerating and spreading relativistic jet 
to broadband data. This model has at most 8 fit parameters: 
Eiso, no, 6q, 9„hs (observer angle), e^., cb, £,n, P- Observer lu- 
minosity distance di and redshift z are assumed to have been 
determined separately. Not all fit parameters need to be in- 
cluded in the fit and any parameter can be fixed to a specific 
value. For example, = 1 and Oohs = rad are commonly 
used. 

The fit code takes as input the full set of broadband data 
points, all expressed in mJy. We use the downhill simplex 
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method dNelder & Meadl[T965l) combined w ith simulated an- 
nealing to minimize dKirkpatrick et al.l ll983; Press et al. 

mWe also use the suggestion from iNelder & MeadI 
to set the result for trial parameters outside of a spec- 
ified parameter domain (e.g. dobs < 0) equal to a very large 
number, which has the effect that the trial will be discarded 
before the downhill simplex iteration has completed. 

The code is written in parallel. The broadband data points 
are distributed over the different computer cores and each core 
calculates the box-based flux counterpart for the data points 
it gets assigned, for each iteration of the fit parameters. Al- 
though the code can also be run on a single core, in practice 
the size of broadband data sets implies that even if the calcula- 
tion of a single datapoint takes mere seconds, the total amount 
of calculation time required for the entire data set can become 
substantial. This is relevant especially in the case of an itera- 
tive fit that requires thousands of iterations (although strongly 
dependent on computer hardware, number of data points and 
numerical accuracy, the procedure can then still take days to 
complete). 

In order to obtain a measure of the error on the fit vari- 
ables, a Monte Carlo (MC) procedure is followed where the 
initial fluxes of the data set are randomly perturbed with an 
amplitude based on their error bars, and with a random num- 
ber drawn from a Gaussian distribution, and then the fit is 
redone. This procedure is repeated a large number of times, 
i.e. 10,000. We take the lowest 68.3% of the resulting val- 
ues and the extremes for the fit parameters within this subset 
determine their Icr uncertainties. 

Although the code allows for an MC calculation where the 
full box calculation is done each time a new model flux is re- 
quired, the amount of time needed for the full MC run can be- 
come prohibitive. In order to circumvent the need for 10,000 
data fits (consisting of thousands of flux calculations per data 
point each), the code offers an alternative approach for esti- 
mating the fit variable errors by calculating the light curves 
for fit variables other than the best fit from a series expansion 
in terms of the fit variables around the best fit. This means that 
instead of a complete flux calculation, at first the best fit result 
is calculated in detail. In addition to this the partial derivatives 
of the flux with respect to the different fit parameters are cal- 
culated. From the base light curve and the derivatives it is 
now possible to estimate the light curves at slightly differing 
values of the fit parameters. These are the values that will in 
practice be probed when is minimized for a perturbed data 
set. 

Rather than the derivative dF/dno (or any fit parameter 
other than p,9obs,0o), we use 91ogF/91ogno for this ap- 
proach. The reason is that we know from analytical model- 
ing that the flux in each spectral regime scales according to 
F oc E^gti'^'ef^eg^i^'^'' etc., where the coefficie nts a, are ei- 
ther c onstant or linearly dependent on p (see e.g. iGranot et al .l 
12001b . The method of calculating logF (rather than F) from 
some base value plus partial derivatives is therefore accurate 
beyond a mere first order approximation (in case of fixed p). 
Otherwise, the accuracy of the series expansion approach is 
set by the deviation of a given set of trial values from the 
best fit values. The maximum for this deviation is ultimately 
determined by the error on the data points. We did not use 
the logarithm of the angles because the observer angle can be 
equal to zero. 

6. GRB 990510: A CASE STUDY 



The "box"-based tool for GRB afterglow modeling pre- 
sented in this paper can be applied to any broadband data set. 
Good spectral and temporal coverage is necessary to accu- 
rately determine all of the macro- and microphysical param- 
eters. The broadband spectrum should span radio to X-ray 
frequencies to encompass all three characteristic frequencies 
of the synchrotron spectrum, and both early and late-time data 
are needed to determine for instance the opening angle of the 
jet and the observer angle. We note that the tool allows for fit- 
ting of limited data sets by fixing some of the fit parameters. 

To demonstrate the capabilities of our tool, we have se- 
lected the afterglow of GRB 990510. There are light curves 
available for this source in X-rays, in various optical bands, 
and at two radio frequencies. Historically, GRB 990510 was 
the first strong case for an achromatic jet break, and the broad- 
band afterglow has been modeled by several authors, all using 
analytical expressions. Here we show the results of our mod- 
eling with RHD simulations. 

The modeling tool requires fluxes in mJy at all frequen- 
cies, which means that some conversions have to be ap- 
plied to the optical and X-ray data. The radio data at 4. 8 
and 8.7 GHz were taken directly from Harrison et all (Il999h . 
The X-ray count rates from Kuulkers et al. (2000) were con- 
verted to mJy using the conversion factors given in that pa- 
per. I n our modeling we ha v e used the V -, R- and 1-ba nd data 
fromlHarrison et all (119991): llsrael et all ( 1999): Stane k et all 
(11999^: IPietrzvnski & Udalskil (ll999l):" lBloom et al.i (11999h : 
Beuermann et al. (1999). We have corrected the optical 
magnitudes for Galactic extinction with E(B -V) = 0.20 
( fSchle gel et al. 1998) before converting the magnitudes into 
fluxes. Starling et al. (2007) have shown that the host galaxy 
extinction is negligible for GRB 990510 based on modeling 
of the X-ray to optical spectrum. 

6.1. Fit results 

We have performed three different fits to the 205 data points 
of the GRB 990510 afterglow, for which we show the re- 
sulting fit parameters in Table |2] and the accompanying best 
fit light curves in Figures [T0l - IT2l Figure [TO] shows the fit 
results for a fit with fixed dots = and = 1, which can 
be d irectly compare d to the broadband fits performed by 
iPanaitescu &lCumarl (|2002|) who use analytical expressions. 
Figure [TT| shows the best fit light curves for a fit with as 
a free parameter and fixed 9obs = 0, and Figure [12] for a fit 
without any fixed parameters. 

Table |2] shows a clear spread in best fit parameters for the 
three fits we performed a nd also some differences with the 
results from Panaitesc u & Kumar! ( 12002 ). These differences 
can be mostly attributed to the fact that the synchrotron self- 
absorption frequency Va is not very well defined by this partic- 
ular data set. The coverage at radio frequencies is fairly sparse 
compared to the optical and X-ray bands, and the flux uncer- 
tainties in the radio are also larger than at higher frequencies. 
For all three of our fits i/g has a value around 10^ Hz at 1 day 
after the burst, while the peak frequency ^ 10'^ Hz at that 
time. The good coverage in the optical and X-ray bands en- 
ables an accurate determination of the cooling frequency Vc, 
which is situated just above the optical bands, and the value of 
p. In contrast to the results from Panaitescu & Kumar (2002), 
we find p > 2 for all three our fits rather than 1.8. Although 
our p > 2 lies within the error bar of their work, the converse 
is not true, which confirms p > 2, and that there is thus no 
need to include an additional high energy cut-off on the rela- 
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TABLE 2 

Best-fit model parameters and 1 a errors. The column labeled PK02 lists the fit results from IPanaitescu & KumarI i I2002I) . 



PK02 



on-axis, fixed on-axis 



off-axis 



00 (rad) 

(erg) 
no (cm"^) 
6obs(rad) 
P 
SB 



5.<i X 10-2 
2.9+j]:|' X 10-1 



xT 



(fixed) 

^•"-'-o.ooe 
5.2+: 



_2 9 X 10 
2.5il X 10-2 
1 (fixed) 



7.5!«-2 X 10-2 
1.8+JJ-3 X 10" 
X 10-2 

(fixed) 
9 70+O.O6 
'^■'^"-O.Ol 

4.6^^1 X 10-3 
3.73^;^« X 10- 

1 (fixed) 



10-2 
10« 
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Fig. 10.— Fit results for the GRB 990510 on-axis fit (6»„6,, = 0) with fixed 
5ai = 1. The reduced x~ for 205 data points and 6 fit parameters is 6.389. 
For clarity of presentation, some fluxes have been multiplied by the indicated 
factors. 



tivistic particle distribution. The value of p has also been de- 
termined by several other authors based on optical and X-ray 
light curve slopes (e.g. Harrison et al. 1999; Kuulk ers et alJ 
HOOO; Panaitescu & Kumar 200]]) 

or a combination of light 
curves and optical to X-ray spectra dStarUng et al.ll200^ . In 
those studies the value for p falls in the range 2.1 -2.2, con- 
sistent with the p value we have obtained in our fit without 
any fixed parameters. 



J-u -1 1 

10 10 10 

observer time (days) 

Fig. 11.— Fit results for GRB 990510 broadband on-axis fit (9ohs = 0). 
The reduced for 205 datapoints and 7 fit parameters is 5.389. For clarity 
of presentation, some fluxes have been multiplied by the indicated factors. 

It is clear from Table |2] that adding extra parameters intro- 
duces quite a strong variation in some of the parameters. In 
particular, the energy and circumburst density change from 
the on-axis to the off-axis fit with more than an order of mag- 
nitude, while in the off-axis fit the opening angle becomes a 
factor of two smaller by introducing a non-zero observer an- 
gle. This shows the importance including the observer angle 
in broadband modeling of GRB afterglows, although this does 
require well-sampled light curves across the whole spectrum. 
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Fig. 12.— Fit results for GRB 990510 broadband off-axis fit, with the 
observer angle included as a fit parameter. The reduced x~ 205 datapoints 
and 8 fit parameters is 3.235. For clarity of presentation, some fluxes have 
been multiplied by the indicated factors. 

More detailed studies of well-sampled broadband afterglows 
will be presented in a future paper. 

7. DISCUSSION 

In this paper we present a method to directly fit light curves 
based on 2D hydrodynamical jet simulations to broadband 
afterglow data. This provides a clear improvement over fits 
based on analytical models, which are not able to take com- 
plex features of the jet dynamics into account, such as the 
radial fluid profile, slow deceleration from ultra-relativistic to 
non-relativistic outflow, the sideways spreading, and resulting 
inhomogeneity along the shock front. The iterative fit proce- 
dure is possible because (/) we have shown that the jet evo- 
lution is scale invariant with respect to explosion energy and 
circumburst medium density, and (//) the results from high- 
resolution parallel RHD simulations can be summarized in a 
very compact form. The compressed version, a 'box' sum- 
mary, of the simulation 'grid' data is possible once the blast 
wave lateral extent, radius and radial width are known from 
the data at each emission time. The predicted flux is calcu- 
lated for each data point for a given set of explosion and ra- 
diation parameters, and because the code takes into account 
both electron cooling and synchrotron self-absorption, it is 
applicable to the entire broadband afterglow spectrum. 



In order to set up the boxes, a total of 19 RHD simulations 
were performed in 2D, with initial opening angles varying be- 
tween 0.045 and 0.5 radians. Light curves for intermediate 
opening angles between different tabulated simulation results 
are obtained via interpolation at the fluid level. The simula- 
tions are run for a long time in order to ensure that late time 
non-relativistic features such as the rise of the counterjet are 
covered as well. At very early times the outflow conforms 
to the self-similar Blandford-McKee solution and this is used 
directly to calculate radiation from emission times before the 
starting point of the simulations. A comparison of the evolu- 
tion of the 19 simulations reveals that, although hard to pre- 
dict analytically, the evolution from strongly collimated BM 
jet to semi-spherical Sedov-Taylor blast wave is smooth. For 
small opening angles, the intermediate stage between the two 
asymptotes is more pronounced. 

We present a number of tests for the resolution of box-based 
light curves by direct comparison to the simulations under- 
lying the box summaries and to earlier work. The earlier 
light curves have been generated directly from simulations 
using different methods: by applying linear radiative trans- 
fer, and by summation of the emitted power (in the optically 
thin case). The differences between box based and simulation 
based light curves are found to be a few percent at most, and 
the box summaries are found to correctly capture the shock 
profile at the different stages of blast wave evolution. 

We have used GRB 990510 as a case study to demonstrate 
our method, because it has been observed over a wide range 
of frequencies. A simultaneous fit has been performed to data 
at two radio frequencies, three optical bands and in X-rays. 
There are some substantial differences between our best fit 
values and earlier fit values found by iPanaitescu & KumaJ 
(120021) . but also between the fits with various parameters 
fixed that we have performed. Most strikingly, including 
a non-zero observing angle in our modeling changes has a 
large impact on the obtained values for the blast wave en- 
ergy and circumburs t densi ty. Furthermore, in contrast with 
IPanaitescu & KuiiimI ( 120021) we find a value for the electron 
energy distribution index p > 2 and we thus do not need to 
include a high energy cut-off on the relativistic particle distri- 
bution. 

The more accurate fits possible by the method of this paper 
help to further constrain the physics that shape the gamma-ray 
burst afterglows. Explosion energy and circumburst density 
provide important clues to the nature of the progenitor star 
and burst environment. In addition, the type of fits to data sets 
covering a long time span that the method makes possible, 
where the complexities of the intermediate stage dynamics 
and the shape of the jet break are fully included, are neces- 
sary to establish a baseline for studies of the effect of more 
detailed models of the microphysics. The evolution of micro- 
physical parameters, su ch as eg, is discussed in various papers 
(Van Eerten et al.ll2010bl: iPanaitescu et al.il2006t iFilgas et aP 
2011). 

Even without utilizing the possibility of fitting broadband 
data directly to simulation results, the fact that the box sum- 
maries cover not just a large number of simulations directly, 
but through scaling and interpolation, exploring the com- 
plete parameter space of impulsive jets in an ISM environ- 
ment (with the inclusion of stellar wind environments being 
straightforward) should prove useful. It will allow us to test 
radiation mechanisms of physical interest other than pure syn- 
chrotron radiation in a realistic context. As long as there 
is no significant feedback on the dynamics of the outflow 
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or dynamically relevant magnetic field involved, any gener- 
alization is possible even if it includes scattering. Exam- 
ples of interest for different radiative processes are given by 
iGiannios & Scitkovskv (2009); PetroDoulou & Mastichiadii 
O009). 

The source code of the broadband fit program 
described in this paper is publicly available at 
http : //cosmo . nyu . edu/ af terglowlibrary 
The code has different settings: not only can it be used to 
fit box based light curves to a data set and establish the 
uncertainties in the best fit parameters, it can also be used 
to generate light curves and spectra directly for arbitrary 
frequencies, observer angles, observer times, and explo- 
sion and radiation parameters. This could be helpful for 
directly exploring how the different regions of the parameter 
space determine the shape of the afterglow, and for quickly 
creating light cur ves that can be expected or looked for in 
surveys (se e e. g. iNakar & PiranI lIOlTt [Roberts et all 120 lit 
iMetzger & Berger 20Tll)! 

A number of practical improvements can be made to the 
code. In the case of very large data sets, even for the paral- 
lel version, where the datapoints to be calculated are evenly 
distributed on the cores, the total calculation time can become 
unwieldy. A remedy to this would be to no longer recalcu- 
late the flux independently for each data point but to calcu- 
late the flux at a fixed (large) number of values for observer 
time and frequency, chosen such to evenly cover the available 
data. The flux at the exact datapoint values can then be de- 



termined by interpolation between these values. Because the 
light curves are smooth, this will not impact the accuracy of 
the fit. For very long data sets (some GRBs have now been ob- 
served for several years, e.g. GRB 030329, Frail et al. 200(1 
Ivan der Horst et al.'2008V more long term simulation data can 
be added to the boxes. Boxes can also be generated for a stel- 
lar wind environment, since the same scale invariances apply. 
These applications and improvements will be presented in fu- 
ture work. 
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APPENDIX 

EMISSION AND ABSORPTION COEFFICIENTS 

The expressions for the emission and absorption coefficient are drawn from lGranot et al.l (11999 1). based in turn on ISari et alj 
119981 and for completeness we provide below the exact forms that we have implemented in our radiation cod43- The code solves 
the radiative transfer equation using 

Al^ = {J^-aMcAt, (Al) 

where At the lab frame time difference between two snapshots, ji, the emission coefficient and the absorption coefficient. 
The code subsequently calculates the flux by integrating over the area A of the emission plane, according to 

F, = -^ jdAI,. (A2) 

Here di is the luminosity distance and z the redshift. Separating the coefficients into frequency dependent and non-frequency 
dependent components, = js x jf and = as x a/, we use for the non-frequency dependent components: 

• ^^"'-S' V3qlip-l)ip + 2) 

Js = 9.6323- 2^7i — S"^' "5 = T7 Tl. ^Nn B -f(l- ^n). (A3) 

Here q,, is the electron charge, the electron mass, n' the comoving number density, B' the comoving magnetic field strength, /3 
the fluid flow velocity as fraction of c, and fi the angle between the outflow and the observer direction. The frequency dependent 
part jf is given by 

if = { W Iv'J'-^'i^ if yL<y' < (A4) 

when i/^ < v'c and by 

r {v'/y'^y" ifiy' <y',<y;„, 

jf={ Wly'cT''^ \fv[<v' <v'„,, (A5) 

I WJy[r"^(y' ly'^yi^ if y[ < yL < y', 

' The terminology used in lVan Eerten et al.l fTOlOaf) is slightly off, where ., m , ■ r u u » j * »u , 

, ' — J T ' ' r-^ , , , J , , explicit as well. The terminology has been corrected tor the cuirent pa- 

he term emissmty has been used to refer to a quantity that should haye been per and now matches -R^^EId a &Ughti£ ^ (M3). The coefficient js in 

labeled emission coefficient, had it not been for the factor 47r that was left | Van Eerten & M acFadven ( 201 ij) was larger by a factor 3/2. 
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otherwise. Here v' denotes the comoving observer frequency. The synchrotron frequency t-,', is given by 



UfB', lL= '-^ ■ (A6) 



Anniec'"" ' \p-l J \£,Nn'meC^ 

The cooling break frequency v'^ is estimated using a global cooling time, leading to 

7^ = ^, t.^t. (A7) 

Note that both and 7^ are comoving with the fluid, while 7^ in lSari et"an (1 1998*) is in the rest frame. 

For the self-absorption coefficient we ignore the effects of cooUng, because in practice the self-absorption break frequency 
^fl ^ ™d because otherwise the limit of accuracy is set in practice by the use of a global cooling time rather than by any 
additional level of detail in the calculation of the absorption coefficient. We have 

SCALE-FREE FLUID EQUATIONS 

In our study we have solved the fluid equations using arbitrary values for explosion energy and density. For reference and 
in order to explictly demonstrate the complete scale-invariance of the simulations we provide the fluid equations in terms of 
dimensionless parameters A, B, 9 below. 

The special relativistic fluid dynamics equations in spherical coordinates, assuming symmetry in angle (f) in the x-y plane 
around the jet axis, are as follows: 

d 8 2 Id 

—lP + (^+ -)lPPr + —^^IPPe sin 61 = 0, 

ct or r r sm u oO 

■^(hj^ -p- pjc^) + ( I- + -)(hj^l3r - p-fc'Pr) + ^ [{hi^fie - Plc^Pe) sin ^] = 0, 
act or r r sm oO 

-hj^P, + (|- + -)(/i7'/3? + P)+ i^h^^PrPe = 0, 
ct Or r r sm oO 

-Wl3e + (|- + -)hj^l3el3r + -^^ihl^Pl + P) = 0- (Bl) 
ct Or r rsmu OO 

Here /3r and /3e are the fluid velocity components in the r and 9 direction respectively in units of c and h the relativistic enthalpy 
density including rest mass energy density. In terms of scale-free parameters A and B, the partial derivatives in r and ct can be 
expressed as 

d _ A d 2B d 
'd^t~~7tdA^~cidB' 
d _l d 5B d 
dr ct dA r dB' 
Combining these with the fluid equations above yields the scale-invariant forms 

. . d ^ d o d 5B d 2n^ 1 d p ^ . ^ ^ 
(-A — + 2B— )7-^ + (t^- — — + -)7-/3.+ T-^^7-/3esm0 = O, 
oA oB po oA A oB A po A sm 9 09 po 

d_ a h^^-p-pjc^ 2 /i7^/3,-p7c^/3, 1 d (hj^/3e- Pic^f3e)sm9 ^ 

^ dA~^ dB' poc^ '^^dA A dB^ a' poc^ A sin 9 d9 pqc^ 

^ dA dB' poc^ dA A dB a' poc^ Asm9 d9 poc^ 

( A^ I 2fi ^ /^'^'^ I r ^ ^ I I 1 ^ h'r'f3l + p _ 

^ dA dB' poc^ ^dA A dB a' pqc^ Asm9 d9 pqc^ ^ ' 

Quantities such as 7p/po etc. are dimensionless and therefore (scale-invariant) functions of the scale-invariant dimensionless 
parameters A, B and 9. For radial outflow and for limiting values of A, the fluid equations further reduce to self-similarity. 
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